//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id$
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#include "MDetectorConstruction.hh"

#include "G4Material.hh"
#include "G4Element.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4OpticalSurface.hh"
#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4Sphere.hh"
#include "G4Colour.hh"

#include <iostream>
#include <fstream>
#include "stdlib.h"
#include <vector>

#include <cmath>
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//
using namespace std;

   G4int MDetectorConstruction::radiator = 3;
   //G4int radiator;
   // 0: C4F10
   // 1: C3F8
   // 2: C2F6
   // 3: CF4
   // 4: Nitrogen
   // 5: Argon
   G4int MDetectorConstruction::window_Material = 2;
   // 0: SiO2
   // 1: MgF2
   // 2: CaF2
   G4double MDetectorConstruction::window_Thickness = 0*mm;	// mm
   G4double N = 2.687e19;	//unit: 1/cm3, the atom density of ideal gas, used to calculate the absorption length
   G4double MDetectorConstruction::H2O_ppm = 0.0001*1e-6;	// should > 0
   G4double MDetectorConstruction::Oxygen_ppm = 0.0001*1e-6;
   G4double PD = 5*mm;	// THGEM + MM = 2*mm, drift = 3*mm
   G4int MDetectorConstruction::MRT_open = 1;	// 0: no optical window, 1: optical window

   
MDetectorConstruction::MDetectorConstruction() 
{
  expHall_x = 0.7*m;
  expHall_y = 0.2*m;
  expHall_z = 0.2*m;
}

MDetectorConstruction::MDetectorConstruction(G4int R, G4int M, G4double T, G4double HP, G4double OP) 
// R: radiator, M: window_Material, T: window_Thickness, HP: H2O_ppm, OP: Oxygen_ppm
{
  expHall_x = 0.7*m;
  expHall_y = 0.2*m;
  expHall_z = 0.2*m;

  MDetectorConstruction::radiator = R;
  MDetectorConstruction::window_Material = M;
  MDetectorConstruction::H2O_ppm = HP*1e-6;	// should > 0
  MDetectorConstruction::Oxygen_ppm = OP*1e-6;
}
MDetectorConstruction::MDetectorConstruction(G4int R, G4int M, G4double T, G4double HP, G4double OP, G4int Mo) 
// R: radiator, M: window_Material, T: window_Thickness, HP: H2O_ppm, OP: Oxygen_ppm
{
  expHall_x = 0.7*m;
  expHall_y = 0.2*m;
  expHall_z = 0.2*m;

  MDetectorConstruction::radiator = R;
  MDetectorConstruction::window_Material = M;
  MDetectorConstruction::H2O_ppm = HP*1e-6;	// should > 0
  MDetectorConstruction::Oxygen_ppm = OP*1e-6;
  MRT_open = Mo;
  cout << "MDetectorConstruction::MRT_open = " << MDetectorConstruction::MRT_open << endl; 
  if (MDetectorConstruction::MRT_open != 1) {MDetectorConstruction::window_Thickness = 0*mm;}
  else { MDetectorConstruction::window_Thickness = T*mm;} 	// mm
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

   
MDetectorConstruction::~MDetectorConstruction(){;}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4VPhysicalVolume* MDetectorConstruction::Construct()
{
//////////////////////////////// Geometry Parameters ////////////////////////////////
   G4double CsI_Thickness = 300*nm;	//nm
   G4double SD_Thickness = PD + window_Thickness + CsI_Thickness;  //mm
   G4double Spherical_Radius = 2*m;
   G4double Focal_Length = Spherical_Radius/2;
   G4double half_Focal_Length = Focal_Length/2;
   G4double Flange_Thickness = 26*mm;
   G4double half_Flange_Thickness = Flange_Thickness/2;
   G4double Mirror_Thickness = 7*mm;
   G4double half_Mirror_Thickness = Mirror_Thickness/2;
   G4double SD_Effective = 20*cm;
   G4double half_SD = SD_Effective/2;
   G4double FR4_Thickness = 8*mm;
   G4double half_FR4_Thickness = FR4_Thickness/2;
   G4double Prototype_Length = Focal_Length + Mirror_Thickness + Flange_Thickness + SD_Thickness + FR4_Thickness - 30/2*mm;
   G4double half_Prototype_Length = Prototype_Length/2;
   G4double Prototype_Width = 22*cm;
   G4double half_Prototype_Width = Prototype_Width/2;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
   //WORLD/////Prototype////**///FR4///**//SD//******************************//Mirror//**//Flange//**////Prototype/////WORLD
             //Prototype////**///FR4///**//SD//******************************//Mirror//**//Flange//**////Prototype//
             //Prototype////**///FR4///**//SD//******************************//Mirror//**//Flange//**////Prototype//
   //WORLD/////Prototype////**///FR4///**//SD//******************************//Mirror//**//Flange//**////Prototype/////WORLD
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


//------------------------------- Materials ------------------------------

  G4int nelements;

  G4double a;  // atomic mass
  G4double z;  // atomic number
  G4double density;
  G4double ncomponents;
  G4double natoms;
//================elements & materials================
  G4Element* C = new G4Element("Carbon", "C", z= 6., a= 12.0107*g/mole);
  G4Element* H = new G4Element("Hydrogen", "H", z= 1., a= 1.00794*g/mole);
  G4Element* O = new G4Element("Oxygen", "O", z= 8., a= 15.9994*g/mole);
  G4Element* F = new G4Element("Fluorine", "F", z= 9., a= 18.998404*g/mole);
  G4Element* Br = new G4Element("Bromine", "Br", z= 35., a= 79.904*g/mole);
  G4Element* Si = new G4Element("Silicon", "Si", z= 14., a= 28.0855*g/mole );
  G4Material* Ar = new G4Material("ArgonGas", z= 18., a= 39.95*g/mole, density= 1.782*mg/cm3);
  G4Material* Al = new G4Material("Aluminium", z= 13., a= 26.982*g/mole, density=26.982*g/cm3);
  G4Material* Fe = new G4Material("Iron", z= 26., a=55.845*g/mole, density= 7.874*g/cm3);
  G4Element* Cs = new G4Element("Cesium", "Cs", z= 55., a=132.90545*g/mole);
  G4Element* I  = new G4Element("Iodine", "I", z= 53., a=126.90447*g/mole);
  G4Element* Mg = new G4Element("Magnesium", "Mg", z= 12., a=24.305*g/mole);
  G4Element* Ca = new G4Element("Calcium", "Ca", z= 20., a=40.078*g/mole);
 
//================FR4================ "epoxy resin"
  G4Material* FR4 = new G4Material("FR4", 1.98281*g/cm3, ncomponents=4);
  FR4->AddElement(C, 38);
  FR4->AddElement(H, 40);
  FR4->AddElement(O, 6);
  FR4->AddElement(Br, 4);
//================C4F10================
  G4double pressure = 1*atmosphere;
  G4double temperature = 293.15*kelvin;
  G4Material* C4F10= new G4Material("C4F10", 11.21*mg/cm3,ncomponents=2, kStateGas, temperature, pressure);
  C4F10->AddElement(C,4);
  C4F10->AddElement(F,10); 
//================C3F8================
  G4Material* C3F8= new G4Material("C3F8", 8.17*mg/cm3,ncomponents=2);
  C3F8->AddElement(C,3);
  C3F8->AddElement(F,8); 
//================C2F6================
  G4Material* C2F6= new G4Material("C2F6", 5.734*mg/cm3,ncomponents=2);
  C2F6->AddElement(C,2);
  C2F6->AddElement(F,6); 
//================CF4================
  G4Material* CF4= new G4Material("CF4", 3.72*mg/cm3,ncomponents=2);
  CF4->AddElement(C,1);
  CF4->AddElement(F,4); 
//================CO2================
  G4Material* CO2 = new G4Material("CO2", 1.977*mg/cm3, ncomponents=2);
  CO2->AddElement(C, 1);
  CO2->AddElement(O, 2);
//================O2================
  G4Material* Oxygen = new G4Material("O2", 1.429*mg/cm3, ncomponents=1);
  Oxygen->AddElement(O, 2);
//================H2O================
  G4Material* H2O = new G4Material("H2O", 1.27*mg/cm3, ncomponents=2);
  H2O->AddElement(H, 2);
  H2O->AddElement(O, 1);
//================Radiator Mixture=================
  G4Material* Radiator_mixture = new G4Material("Radiator_mixture", 6.1*mg/cm3, ncomponents = 3);
  G4Material* Radiator_base;
  switch (radiator)
  {
    case 0: Radiator_base = C4F10; break;
    case 1: Radiator_base = C3F8; break;
    case 2: Radiator_base = C2F6; break;
    case 3: Radiator_base = CF4; break;
  }
  Radiator_mixture -> AddMaterial(Radiator_base, (100 - H2O_ppm*1e-6 - Oxygen_ppm*1e-6)*perCent);
  Radiator_mixture -> AddMaterial(H2O, H2O_ppm*1e-6*perCent);
  Radiator_mixture -> AddMaterial(Oxygen, Oxygen_ppm*1e-6*perCent);

//================Avalanche gas(Ar + CO2)================
  G4Material* Avalanche_gas = new G4Material("Avalanche_gas", 1.79565*mg/cm3, ncomponents = 2);
  Avalanche_gas->AddMaterial(CO2, 7*perCent);
  Avalanche_gas->AddMaterial(Ar, 93*perCent);
//================SiO2===============
  G4Material* SiO2 = new G4Material("SiO2", 2.65*g/cm3,ncomponents=2);
  SiO2->AddElement(Si, 1);
  SiO2->AddElement(O, 2);
//================Quartz================
  G4Material* Quartz = new G4Material("Quartz", 2.649*g/cm3,ncomponents=2);
  Quartz->AddElement(Si, 1);
  Quartz->AddElement(O, 2);
//================MgF2================
  G4Material* MgF2 = new G4Material("MgF2", 3.148*g/cm3, nelements=2);
  MgF2->AddElement(Mg, 1);
  MgF2->AddElement(F, 2);
//================CaF2================
  G4Material* CaF2 = new G4Material("CaF2", 3.8*g/cm3, nelements=2);
  CaF2->AddElement(Ca, 1);
  CaF2->AddElement(F, 2);
//================CsI================
  G4Material* CsI = new G4Material("CsI", 4.53*g/cm3, nelements=2);
  CsI->AddElement(Cs, 1);
  CsI->AddElement(I, 1);

//------------------------------ Generate & Add Material Properties Table ------------------------------
  const G4int nEntries = 45;
//================optical photon energy=================
    G4double PhotonEnergy[nEntries] = 
	{ 2.96171*eV, 3.12387*eV, 3.28378*eV, 3.44595*eV,	
	  3.60586*eV, 3.76577*eV, 3.92793*eV, 4.09009*eV,	
	  4.25225*eV, 4.41441*eV, 
	  4.57432*eV, 4.73423*eV, 4.8964 *eV, 5.05631*eV,	
	  5.21847*eV, 5.37838*eV, 5.53829*eV, 5.6982 *eV, 
	  5.85586*eV, 6.01351*eV,	
	  6.17568*eV, 6.33559*eV, 6.4955 *eV, 6.65541*eV,	
	  6.81532*eV, 6.97523*eV, 7.13514*eV, 7.29505*eV,	
	  7.45721*eV, 7.61486*eV,	
	  7.77252*eV, 7.93243*eV, 8.08784*eV, 8.24324*eV,	
	  8.4009 *eV, 8.55405*eV, 8.71396*eV, 8.86937*eV,	
	  9.01802*eV, 9.17568*eV,	
	  9.32883*eV, 9.48649*eV, 9.63964*eV, 9.79054*eV,	
	  9.9527 *eV 
	}; // 3eV to 10eV, 120nm to 400nm
    G4double PEnergy[nEntries] = 
	{ 2.96171, 3.12387, 3.28378, 3.44595,	
	  3.60586, 3.76577, 3.92793, 4.09009,	
	  4.25225, 4.41441, 
	  4.57432, 4.73423, 4.8964 , 5.05631,	
	  5.21847, 5.37838, 5.53829, 5.6982 , 
	  5.85586, 6.01351,	
	  6.17568, 6.33559, 6.4955 , 6.65541,	
	  6.81532, 6.97523, 7.13514, 7.29505,	
	  7.45721, 7.61486,	
	  7.77252, 7.93243, 8.08784, 8.24324,	
	  8.4009 , 8.55405, 8.71396, 8.86937,	
	  9.01802, 9.17568,	
	  9.32883, 9.48649, 9.63964, 9.79054,	
	  9.9527  
	}; // without unit

//))))))))))))))))Radiator((((((((((((((((((((((((((	
//================C4F10 refractive index================
	G4double RefractiveIndexC4F10[nEntries] = 
	{ 1.00137295, 1.00137585, 1.00137971, 1.00138454, 
	  1.0013913,  1.0013971,  1.0014,     1.00140097,
	  1.00140386, 1.0014087,  
	  1.00141739, 1.00142705, 1.00143478, 1.00143961,
	  1.00144541, 1.00145507, 1.0014657,  1.00147246,
	  1.00148213, 1.00149179,
	  1.00149758, 1.00150435, 1.00151304, 1.00152464,
	  1.0015372,  1.00154879, 1.00156039, 1.00157005,
	  1.00157778, 1.0015913,  
	  1.00160483, 1.00161256, 1.00163188, 1.00164928,
	  1.00166473, 1.00168309, 1.00169082, 1.00170821,
	  1.00173527, 1.001743,   
	  1.00176135, 1.00178068, 1.00180386, 1.00182512,
	  1.00182609
	};

//================C4F10 physics table==================
  G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
  myMPT1->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexC4F10,nEntries)
        ->SetSpline(true);
  C4F10->SetMaterialPropertiesTable(myMPT1);

//================C3F8 refractive index================
 G4double RefractiveIndexC3F8[nEntries];
 for (G4int i = 0; i < nEntries; i++)
 {
   G4double wavelength = 1240/PEnergy[i];
   G4double refractive = 1 + 0.279337e-6/(2.408623e-4 - pow(wavelength, -2));
   RefractiveIndexC3F8[i] = refractive;
 }

//================C3F8 physics table==================
  G4MaterialPropertiesTable* myMPT_C3F8 = new G4MaterialPropertiesTable();
  myMPT_C3F8->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexC3F8,nEntries)
        ->SetSpline(true);
  C3F8->SetMaterialPropertiesTable(myMPT_C3F8);

//================C2F6 refractive index================
  G4double RefractiveIndexC2F6[nEntries];
  for (G4int i = 0; i < nEntries; i++)
  {
    G4double wavelength = 1240/PEnergy[i];
    G4double refractive = 1 + 0.18994e-6/(2.33313e-4 - pow(wavelength, -2));
    RefractiveIndexC2F6[i] = refractive;
  }

//================C2F6 physics table==================
  G4MaterialPropertiesTable* myMPT_C2F6 = new G4MaterialPropertiesTable();
  myMPT_C2F6->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexC2F6,nEntries)
        ->SetSpline(true);
  C2F6->SetMaterialPropertiesTable(myMPT_C2F6);

//================CF4 refractive index================
  G4double RefractiveIndexCF4[nEntries];
  for (G4int i = 0; i < nEntries; i++)
  {
    G4double wavelength = 1240/PEnergy[i];
    G4double refractive = 1 + (0.124523e-6/(2.61154e-4 - pow(wavelength, -2)))*106/101;
    RefractiveIndexCF4[i] = refractive;
  }

//================CF4 physics table==================
  G4MaterialPropertiesTable* myMPT_CF4 = new G4MaterialPropertiesTable();
  myMPT_CF4->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexCF4,nEntries);
        //->SetSpline(true);
  CF4->SetMaterialPropertiesTable(myMPT_CF4);

//================H2O physics table=================
  G4int num_H2O = 0;
  fstream H2O_absorption;
  H2O_absorption.open("../Reference/Absorption/water_absorption_PHENIX.csv", ios::in);
  if (H2O_absorption.is_open()) {G4cout << "H2O_absorption file opened!" << G4endl;}
  else {G4cout << "H2O_absorption open failed! " << G4endl;exit(1);}  

  vector<double> H2O_wavelength, H2O_AbsorptionIndex;
  
  while (!H2O_absorption.eof())
  {
    G4double buffer1 = 0, buffer2 = 0;
    H2O_absorption >> buffer1 >> buffer2;
    H2O_wavelength.push_back(buffer1);
    H2O_AbsorptionIndex.push_back(buffer2);
    num_H2O++;
  }
  num_H2O = num_H2O - 1;
  G4double H2O_Ephoton[num_H2O];
  G4double H2O_Absorption[num_H2O]; 
  for (int i = 0; i < num_H2O; i++)
  {
    H2O_Ephoton[i] = 1240/H2O_wavelength[i]*eV;
    H2O_Absorption[i] = 1/(H2O_AbsorptionIndex[i]*N*H2O_ppm*1e-18)*cm;	//unit: cm, 1e-6 is from ppm, 1e-18 is from cross section
  }

//================H2O physics table================
  G4MaterialPropertiesTable* myMPT_H2O = new G4MaterialPropertiesTable();
  myMPT_H2O->AddProperty("ABSLENGTH", H2O_Ephoton, H2O_Absorption,num_H2O);
        //->SetSpline(true);
  H2O->SetMaterialPropertiesTable(myMPT_H2O);

//================Oxygen physics table=================
  G4int num_Oxygen = 0;
  fstream Oxygen_absorption;
  Oxygen_absorption.open("../Reference/Absorption/oxygen_absorption_PHENIX.csv", ios::in);
  if (Oxygen_absorption.is_open()) {G4cout << "Oxygen_absorption file opened!" << G4endl;}
  else {G4cout << "Oxygen_absorption open failed! " << G4endl;exit(1);}  

  vector<double> Oxygen_wavelength, Oxygen_AbsorptionIndex;
  
  while (!Oxygen_absorption.eof())
  {
    G4double buffer1 = 0, buffer2 = 0;
    Oxygen_absorption >> buffer1 >> buffer2;
    Oxygen_wavelength.push_back(buffer1);
    Oxygen_AbsorptionIndex.push_back(buffer2);
    num_Oxygen++;
  }
  num_Oxygen = num_Oxygen - 1;
  G4double Oxygen_Ephoton[num_Oxygen];
  G4double Oxygen_Absorption[num_Oxygen]; 
  for (int i = 0; i < num_Oxygen; i++)
  {
    Oxygen_Ephoton[i] = 1240/Oxygen_wavelength[i]*eV;
    Oxygen_Absorption[i] = 1/(Oxygen_AbsorptionIndex[i]*N*Oxygen_ppm*1e-18)*cm;	//unit: cm, 1e-6 is from ppm, 1e-18 is from cross section
  }

//================Oxygen physics table================
  G4MaterialPropertiesTable* myMPT_Oxygen = new G4MaterialPropertiesTable();
  myMPT_Oxygen->AddProperty("ABSLENGTH", Oxygen_Ephoton, Oxygen_Absorption,num_Oxygen);
        //->SetSpline(true);
  Oxygen->SetMaterialPropertiesTable(myMPT_Oxygen);

//+++++++++++++++++sum of all contaminants++++++++++++
  G4double Total_Ephoton[num_Oxygen], Total_Absorption[num_Oxygen];	// in case num_Oxygen = num_H2O
  for (int i = 0; i < num_Oxygen; i++)
  {
    Total_Ephoton[i] = Oxygen_Ephoton[i];
    //Total_Absorption[i] = Oxygen_Absorption[i] + H2O_Absorption[i];
    Total_Absorption[i] = 1/((Oxygen_AbsorptionIndex[i]*Oxygen_ppm + H2O_AbsorptionIndex[i]*H2O_ppm)*N*1e-18);
  }

//================Radiator physics table==============
  G4MaterialPropertiesTable* myMPT_Radiator_mixture = new G4MaterialPropertiesTable();
  myMPT_Radiator_mixture->AddProperty("ABSLENGTH", Total_Ephoton, Total_Absorption,num_Oxygen);
        //->SetSpline(true);
  switch (radiator)
  {
    case 0: myMPT_Radiator_mixture->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexC4F10,nEntries)->SetSpline(true); break;
    case 1: myMPT_Radiator_mixture->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexC3F8,nEntries)->SetSpline(true); break;
    case 2: myMPT_Radiator_mixture->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexC2F6,nEntries)->SetSpline(true); break;
    case 3: myMPT_Radiator_mixture->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexCF4,nEntries);break;//->SetSpline(true); break;
  }
  Radiator_mixture->SetMaterialPropertiesTable(myMPT_Radiator_mixture);
 

//))))))))))))))))Optical Window(((((((((((((((((((
//================SiO2================
  //G4cout << "^^^^^^^^^^^^^^^^^^^^SiO2^^^^^^^^^^^^^^^^^^^^" << G4endl;
  G4double RefractiveIndexSiO2[nEntries], PhotonEnergy_SiO2[nEntries];
  for (G4int i = 0; i < nEntries; i++)
  {
    G4double wavelength = 1.240/PEnergy[i];
    G4double refractive = sqrt(1 + 0.6961663*pow(wavelength,2)/(pow(wavelength,2) - pow(0.0684043,2)) + 0.4079426*pow(wavelength,2)/(pow(wavelength,2)-pow(0.1162414,2)) + 0.8974794*pow(wavelength,2)/(pow(wavelength,2) - pow(9.896161,2)));
    RefractiveIndexSiO2[i] = refractive;
  }

//================SiO2 physics table================
  G4MaterialPropertiesTable* myMPT_SiO2 = new G4MaterialPropertiesTable();
  myMPT_SiO2->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexSiO2,nEntries);
        //->SetSpline(true);
  SiO2->SetMaterialPropertiesTable(myMPT_SiO2);

//================MgF2================
  //G4cout << "^^^^^^^^^^^^^^^^^^^^MgF2^^^^^^^^^^^^^^^^^^^^" << G4endl;
  G4double RefractiveIndexMgF2[nEntries], PhotonEnergy_MgF2[nEntries];
  for (G4int i = 0; i < nEntries; i++)
  {
    G4double wavelength = 1.240/PEnergy[i];
    G4double refractive = sqrt(1 + 0.27620 + 0.60967*pow(wavelength,2)/(pow(wavelength,2) - pow(0.08636,2)) + 0.0080*pow(wavelength,2)/(pow(wavelength,2)-pow(18.0,2)) + 2.14973*pow(wavelength,2)/(pow(wavelength,2) - pow(25.0,2)));
    RefractiveIndexMgF2[i] = refractive;
  }

//================MgF2 physics table================
  G4MaterialPropertiesTable* myMPT_MgF2 = new G4MaterialPropertiesTable();
  myMPT_MgF2->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexMgF2,nEntries);
        //->SetSpline(true);
  MgF2->SetMaterialPropertiesTable(myMPT_MgF2);

//================CaF2================
  //G4cout << "^^^^^^^^^^^^^^^^^^^^CaF2^^^^^^^^^^^^^^^^^^^^" << G4endl;
  G4double RefractiveIndexCaF2[nEntries], PhotonEnergy_CaF2[nEntries];
  for (G4int i = 0; i < nEntries; i++)
  {
    G4double wavelength = 1.240/PEnergy[i];
    G4double refractive = sqrt(1 + 0.443749998*pow(wavelength,2)/(pow(wavelength,2) - 0.00178027854) + 0.444930066*pow(wavelength,2)/(pow(wavelength,2)- 0.00788536061) + 0.150133991*pow(wavelength,2)/(pow(wavelength,2) - 0.0124119491) + 8.85319946*pow(wavelength,2)/(pow(wavelength,2) - 2752.28175));
    RefractiveIndexCaF2[i] = refractive;
  }

//================CaF2 physics table================
  G4MaterialPropertiesTable* myMPT_CaF2 = new G4MaterialPropertiesTable();
  myMPT_CaF2->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexCaF2,nEntries);
        //->SetSpline(true);
  CaF2->SetMaterialPropertiesTable(myMPT_CaF2);


//================CsI================
  G4double RefractiveIndex3[nEntries] = 
	    {2.30,  2.30,  2.30,  2.25,
	     2.25,  2.21,  2.14,  2.07,
	     2.04,  2.01,
	     1.98,  1.96,  1.94,  1.92,
	     1.91,  1.89,  1.882, 1.87,
	     1.863, 1.858,
             1.852, 1.845, 1.839, 1.834, 
	     1.830, 1.823, 1.820, 1.818, 
	     1.812, 1.810, 
   	     1.806, 1.804, 1.801, 1.799,
  	     1.796, 1.794, 1.793, 1.790,
	     1.789, 1.787,
	     1.785, 1.783, 1.781, 1.679,
	     1.677};

  G4double Absorption3[nEntries] = 
	     {0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm, 
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm, 
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
              0.01*nm, 0.01*nm,
	      0.01*nm, 0.01*nm, 0.01*nm, 0.01*nm,
	      0.01*nm};
//================CsI physics table================
  G4MaterialPropertiesTable* myMPT3 = new G4MaterialPropertiesTable();
  myMPT3->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex3,nEntries);
        //->SetSpline(true);
  myMPT3->AddProperty("ABSLENGTH", PhotonEnergy, Absorption3,nEntries);
        //->SetSpline(true);
  CsI->SetMaterialPropertiesTable(myMPT3);

//================Avalanche gas================
  G4double RefractiveIndex4[nEntries] =
            {1.,  1.,  1.,  1.,
	     1.,  1.,  1.,  1.,
	     1.,  1., 
	     1.,  1.,  1.,  1.,
	     1.,  1.,  1.,  1.,
  	     1.,  1.,
	     1.,  1.,  1.,  1.,
	     1.,  1.,  1.,  1., 
	     1.,  1.,
	     1.,  1.,  1.,  1.,
	     1.,  1.,  1.,  1.,
             1.,  1.,
	     1.,  1.,  1.,  1., 
	     1.};
  G4double Absorption4[nEntries] = 
	     {1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm, 
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm, 
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
              1E10*nm, 1E10*nm,
	      1E10*nm, 1E10*nm, 1E10*nm, 1E10*nm,
	      1E10*nm};

//================Avalanche gas physics table================
  G4MaterialPropertiesTable* myMPT4 = new G4MaterialPropertiesTable();
 // myMPT4->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex4,nEntries)
 //       ->SetSpline(true);
  myMPT4->AddProperty("ABSLENGTH", PhotonEnergy, Absorption4,nEntries);
        ///->SetSpline(true);
  Avalanche_gas->SetMaterialPropertiesTable(myMPT4);

// ------------------------------- Volumes --------------------------------
//================= The experimental Hall================
  G4Box* expHall_box = new G4Box("World",expHall_x,expHall_y,expHall_z);
  G4LogicalVolume* expHall_log;
  expHall_log = new G4LogicalVolume(expHall_box, Avalanche_gas,"World",0,0,0); 
  //expHall_log = new G4LogicalVolume(expHall_box, C3F8,"World",0,0,0);  break;
  //expHall_log = new G4LogicalVolume(expHall_box, C2F6,"World",0,0,0);  break;
  //expHall_log = new G4LogicalVolume(expHall_box, CF4,"World",0,0,0);   break;
  //expHall_log = new G4LogicalVolume(expHall_box, Radiator_mixture,"World",0,0,0);
  G4VPhysicalVolume* expHall_phys = new G4PVPlacement(0,G4ThreeVector(),expHall_log,"World",0,false,0);

  // G4Box::::Construct a box with name, and half lengths pX,pY,pZ 
  // G4VPhysicalVolume::::The center of the object is then placed at tlate in the new coordinates of mother volume. 


//================Prototype Section================
  G4Box* Prototype = new G4Box("Prototype", half_Prototype_Length, half_Prototype_Width, half_Prototype_Width);
  G4LogicalVolume* Prototype_log;
  switch (radiator)
  {
    case 0:  Prototype_log = new G4LogicalVolume(Prototype, C4F10,"Prototype_log",0,0,0); break;
    case 1:  Prototype_log = new G4LogicalVolume(Prototype, C3F8,"Prototype_log",0,0,0);  break;
    case 2:  Prototype_log = new G4LogicalVolume(Prototype, C2F6,"Prototype_log",0,0,0);  break;
    case 3:  Prototype_log = new G4LogicalVolume(Prototype, CF4,"Prototype_log",0,0,0);   break;
  }
  Prototype_log = new G4LogicalVolume(Prototype, Radiator_mixture,"Prototype_log",0,0,0); 

  G4double Prototype_Pos_x = 0.0*m;
  G4double Prototype_Pos_y = 0.0*m;
  G4double Prototype_Pos_z = 0.0*m;
  G4VPhysicalVolume* Prototype_phys = new G4PVPlacement(0,G4ThreeVector(Prototype_Pos_x,Prototype_Pos_y,Prototype_Pos_z),Prototype_log,"Prototype_phys",expHall_log,false,0);

//================FR4===============  
  G4Box* SD_FR4 = new G4Box("SD_FR4", half_FR4_Thickness, half_SD, half_SD);
  G4LogicalVolume* FR4_log = new G4LogicalVolume(SD_FR4, FR4, "FR4_log",0,0,0);
  G4double FR4_Pos_x = -1*half_Focal_Length - SD_Thickness - half_FR4_Thickness;
  G4double FR4_Pos_y = 0.0*m;
  G4double FR4_Pos_z = 0.0*m;
  G4VPhysicalVolume* FR4_phys = new G4PVPlacement(0,G4ThreeVector(FR4_Pos_x,FR4_Pos_y,FR4_Pos_z),FR4_log,"FR4_phys",Prototype_log,false,0);

//================SD(sensitive detector)================
  G4Box* SD = new G4Box("Sensitive Detector", SD_Thickness/2, half_SD, half_SD);
  G4LogicalVolume* SD_log = new G4LogicalVolume(SD, Avalanche_gas, "SD_log");
  G4double SD_Pos_x = -1*half_Focal_Length - SD_Thickness/2;	//m
  G4double SD_Pos_y = 0*cm;
  G4double SD_Pos_z = 0*cm;
  G4VPhysicalVolume* SD_phys = new G4PVPlacement(0, G4ThreeVector(SD_Pos_x, SD_Pos_y, SD_Pos_z), SD_log, "SD_phys", Prototype_log, false, 0);

//================MPGD================
  G4Box* MPGD = new G4Box("SD_MicroMegas", 1*mm, half_SD, half_SD);	// for the thickness of MPGD = 2*mm
  G4LogicalVolume* MPGD_log = new G4LogicalVolume(MPGD,FR4,"MPGD_log" );
  G4double MPGD_Pos_x = -SD_Thickness/2 + 1.1*mm;	// for the thickness of MPGD = 2*mm
  G4double MPGD_Pos_y = 0.*m;
  G4double MPGD_Pos_z = 0.*m;
  G4VPhysicalVolume* MPGD_phys = new G4PVPlacement(0, G4ThreeVector(MPGD_Pos_x,MPGD_Pos_y,MPGD_Pos_z), MPGD_log, "MPGD_phys", SD_log, false,0);

//================CsI--photon-electron-converter================
  G4cout << "Window_Thickness = " << window_Thickness << G4endl;
  G4Box* PEC;
  G4LogicalVolume* PEC_log;
  G4VPhysicalVolume* PEC_phys;
  //if (MRT_open == 1)
  //{
    PEC = new G4Box("PEC", CsI_Thickness/2, half_SD, half_SD);	// 300nm thick
    PEC_log = new G4LogicalVolume(PEC, CsI, "PEC_log",0,0,0);
    G4double PEC_Pos_x = SD_Thickness/2 - window_Thickness - CsI_Thickness/2;
    G4double PEC_Pos_y = 0;
    G4double PEC_Pos_z = 0;
    PEC_phys = new G4PVPlacement(0,G4ThreeVector(PEC_Pos_x,PEC_Pos_y,PEC_Pos_z),PEC_log,"PEC_phys",SD_log,false,0);
  //}

//================Optical Window & transmission================
//==========================Optical Window=====================
  //G4cout << "MRT = " << window_Thickness << G4endl;
  G4Box * MRT;
  G4LogicalVolume* MRT_log;
  G4VPhysicalVolume* MRT_phys;
  if (MRT_open == 1)
  {
    MRT = new G4Box("MRT", window_Thickness/2, half_SD, half_SD);
    switch(window_Material)
    {
      case 0: MRT_log = new G4LogicalVolume(MRT, SiO2, "MRT_Log",0,0,0); break;
      case 1: MRT_log = new G4LogicalVolume(MRT, MgF2, "MRT_Log",0,0,0); break;
      case 2: MRT_log = new G4LogicalVolume(MRT, CaF2, "MRT_Log",0,0,0); break;
    }
    G4double MRT_Pos_x = SD_Thickness/2 - window_Thickness/2;
    G4double MRT_Pos_y = 0;
    G4double MRT_Pos_z = 0;
    MRT_phys = new G4PVPlacement(0,G4ThreeVector(MRT_Pos_x,MRT_Pos_y,MRT_Pos_z),MRT_log,"MRT_log",SD_log,false,0);
  }


//================Spherical Reflection Film================
  G4Sphere* Spherical_Film = new G4Sphere("Spherical-Film", Spherical_Radius-100*nm, Spherical_Radius, 0, 2*3.1415, 0, 3.1415);//-0.05, 0.10, 1.5199, 0.10);
  G4LogicalVolume* SF_log = new G4LogicalVolume(Spherical_Film, Al, "SF_log", 0,0,0);
  G4double SF_Pos_x = -3*half_Focal_Length;
  G4double SF_Pos_y = 0*m;
  G4double SF_Pos_z = 0*m;
  G4RotationMatrix* fieldRot = new G4RotationMatrix();
  fieldRot->rotateY(90.*deg);
  G4VPhysicalVolume* SF_phys = new G4PVPlacement(fieldRot, G4ThreeVector(SF_Pos_x, SF_Pos_y, SF_Pos_z), SF_log, "SF_phys", Prototype_log, false, 0);

//================Spherical Shell Section================
  G4Sphere* Spherical_Mirror = new G4Sphere("Spherical-Mirror", Spherical_Radius, Spherical_Radius+Mirror_Thickness, -0.05, 0.10, 1.5199, 0.10);
  //G4Sphere* Spherical_Mirror = new G4Sphere("Spherical-Mirror", 2*m, 2.003*m, -0.03, 0.06, 1.540796, 0.06);
  //G4Sphere* Spherical_Mirror = new G4Sphere("Spherical-Mirror", 2*m, 2.003*m, -3.1415926/2, 3.1415926, 0, 3.1415926);
  G4LogicalVolume* SM_log = new G4LogicalVolume(Spherical_Mirror, Quartz, "SM_log", 0,0,0);
  G4double SM_Pos_x = -3*half_Focal_Length;
  G4double SM_Pos_y = 0*m;
  G4double SM_Pos_z = 0*m;
  G4VPhysicalVolume* SM_phys = new G4PVPlacement(0, G4ThreeVector(SM_Pos_x, SM_Pos_y, SM_Pos_z), SM_log, "SM_phys", Prototype_log, false, 0);


//================Flange Section================
  G4Sphere* Flange = new G4Sphere("Flange", Spherical_Radius+Mirror_Thickness, Spherical_Radius+Mirror_Thickness+Flange_Thickness,  -0.05, 0.10, 1.5199, 0.10);
  G4LogicalVolume* Flange_log = new G4LogicalVolume(Flange, Fe, "Flange_log",0,0,0);
  G4double Flange_Pos_x = -3*half_Focal_Length;
  G4double Flange_Pos_y = 0.0*m;
  G4double Flange_Pos_z = 0.0*m;
  G4VPhysicalVolume* Flange_phys = new G4PVPlacement(0,G4ThreeVector(Flange_Pos_x,Flange_Pos_y,Flange_Pos_z),Flange_log,"Flange_phys",Prototype_log,false,0);


// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Surfaces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

//================Spherical Shell/Radiator Border Surface================ 
  G4OpticalSurface* OpSMSurface = new G4OpticalSurface("SMSurface");
  OpSMSurface->SetType(dielectric_metal);
  OpSMSurface->SetFinish(polished);
  OpSMSurface->SetModel(glisur);
  G4LogicalSkinSurface* SMSurface = new G4LogicalSkinSurface("SMSurface", SF_log, OpSMSurface);
  G4OpticalSurface* opticalSurface1 = dynamic_cast <G4OpticalSurface*>
        (SMSurface->GetSurface(SF_log)->GetSurfaceProperty());

  if (opticalSurface1) opticalSurface1->DumpInfo();            

  return expHall_phys;
//================Radiator/Oprtical Window Border Surface================ 
  G4OpticalSurface* OpCMRTSurface;
  G4OpticalSurface* OpMCSurface;
  if (MRT_open == 1)
  {
    OpCMRTSurface = new G4OpticalSurface("CMRTSurface");
    OpCMRTSurface->SetType(dielectric_dielectric);
    OpCMRTSurface->SetFinish(polished);
    OpCMRTSurface->SetModel(glisur);
    G4LogicalSkinSurface* CMRTSurface = new G4LogicalSkinSurface("CMRTSurface", MRT_log, OpCMRTSurface);
    G4OpticalSurface* opticalSurface2 = dynamic_cast <G4OpticalSurface*>
          (CMRTSurface->GetSurface(MRT_log)->GetSurfaceProperty());

    if (opticalSurface2) opticalSurface2->DumpInfo();
//================CsI/Optical window Border Surface================ 
  OpMCSurface = new G4OpticalSurface("MCSurface");
  OpMCSurface->SetType(dielectric_metal);
  OpMCSurface->SetFinish(polished);
  OpMCSurface->SetModel(glisur);
  G4LogicalSkinSurface* MCSurface = new G4LogicalSkinSurface("MCSurface", PEC_log, OpMCSurface);
  G4OpticalSurface* opticalSurface3 = dynamic_cast <G4OpticalSurface*>
        (MCSurface->GetSurface(PEC_log)->GetSurfaceProperty());

  if (opticalSurface3) opticalSurface3->DumpInfo();
  }
//**************************************************************************
// Generate & Add Material Properties Table attached to the optical surfaces
//================ Radiator/Spherical Mirror Surface================
  G4int numb = 0;
  //G4cout << "^^^^^^^^^^^^^^^^^^^^Al VUV Reflection^^^^^^^^^^^^^^^^^^^^" << G4endl;
  fstream Al_VUV_reflection;
  Al_VUV_reflection.open("../Reference/Reflectivity/Al-VUV-reflection.csv", ios::in);
  if (Al_VUV_reflection.is_open()) {G4cout << "Al_VUV_reflection file opened!" << G4endl;}
  else {G4cout << "Al_VUV_reflection open failed! " << G4endl;exit(0);}  

  vector<double> wavelength, Al_VUV_ReflectionIndex;
  G4double Reflsum = 0;
  
  while (!Al_VUV_reflection.eof())
  {
    G4double buffer1 = 0, buffer2 = 0;
    Al_VUV_reflection >> buffer1 >> buffer2;
    wavelength.push_back(buffer1);
    Al_VUV_ReflectionIndex.push_back(buffer2);
    Reflsum += buffer2;
    numb++;
  }
  const G4int num = numb - 1;
  G4double Ephoton[num];
  G4double Reflectivity[num]; 
  for (int i = 0; i < num; i++)
  {
    Ephoton[i] = 1.240/wavelength[i]*eV;
    Reflectivity[i] = Al_VUV_ReflectionIndex[i];
    //Reflectivity[i] = 1;
    //Reflectivity[i] = Reflsum/num;
    //cout << "Spherical Mirror Reflectivity[" << i << "] = " << Reflectivity[i] << endl;
  }

  G4double Efficiency[2]   = {0.8, 1.0};
  //G4double Efficiency[num]   = {0.0, 0.73};
  //G4double Efficiency[num] = {1.0, 1.0};

  G4MaterialPropertiesTable *myST1 = new G4MaterialPropertiesTable();

  myST1->AddProperty("REFLECTIVITY", Ephoton, Reflectivity, num);
  //myST1->AddProperty("EFFICIENCY",   Ephoton, Efficiency, num);
  // &&&&&&&&&&&&&&&& "EFFICIENCY" valid only to dielectric_metal surfaces &&&&&&&&&&&&&&&&&&&&&&&

  OpSMSurface->SetMaterialPropertiesTable(myST1);

//================Radiator/Optical Window Surface================
//  G4double Refraction1[num] = 
//		{1.2, 1.2};
//
//  G4double Efficiency1[num] = 
//		{1.0, 1.0};
//
//
//  G4MaterialPropertiesTable *myST2 = new G4MaterialPropertiesTable();
//
//  myST2->AddProperty("RINDEX", Ephoton, Refraction1, num);
//  myST2->AddProperty("EFFICIENCY",   Ephoton, Efficiency1,   num);
//
//  OpCMRTSurface->SetMaterialPropertiesTable(myST2);
//
//================Optical Window/CsI Surface================
  G4MaterialPropertiesTable *myST3;
  if (MRT_open == 1)
  {
    G4double Refraction2[2] = 
          	{1.2, 1.2};
    G4double Efficiency2[2]  =
          	{1.0, 1.0};
    myST3 = new G4MaterialPropertiesTable();

    myST3->AddProperty("RINDEX", Ephoton, Refraction2, num);
    myST3->AddProperty("EFFICIENCY",   Ephoton, Efficiency2,   num);
    OpMCSurface->SetMaterialPropertiesTable(myST3);
  }

  Al_VUV_reflection.close();
//always return the physical World
  return expHall_phys;
}

G4int MDetectorConstruction::GetRadiatorIndex()
{
  return radiator;
}

G4int MDetectorConstruction::GetWindowMaterialIndex()
{
  return window_Material;
}

G4double MDetectorConstruction::GetWindowThickness()
{
  return window_Thickness;
}
